Análise Multivariada

Lab: Análise de Componentes Principais em R

Prof. Washington Santos da Silva

IFMG - Campus Formiga

19 de outubro de 2023

Ambiente Virtual de Aprendizagem

  • Lembre-se que a sala virtual da disciplina está disponível no AVA (Moodle) do IFMG:

  • E que para acessar a sala, basta utilizar as mesmas informações (usuario/senha) que vocês utilizam para acessar o sistema acadêmico Meu IFMG

Summário: Aula em 19/10

R Notebook

  • Por favor criem um R Notebook no RStudio e o nomeiem como pca.Rmd

  • Você pode criar um novo notebook no RStudio usando o menu principal:

    • File -> New File -> R Notebook

  • Trechos de código podem ser inseridos usando o seguinte atalho do teclado:

    • Ctrl + Alt + I (Windows/Linux)
    • or pelo menu Insert a new code chunk.

  • Vamos criar juntos o R Notebook para esta aula.

Dados: USArrests

  • Podemos visualizar a documentação dos dados executando o seguinte comando no console:
?USArrests
  • Analisando a estrutura dos dados:
str(USArrests)
'data.frame':   50 obs. of  4 variables:
 $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
 $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
 $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
 $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...

Dados: USArrests

  • As linhas da data frame contêm os 50 estados, em ordem alfabética:
estados <- row.names(USArrests)
estados
 [1] "Alabama"        "Alaska"         "Arizona"        "Arkansas"      
 [5] "California"     "Colorado"       "Connecticut"    "Delaware"      
 [9] "Florida"        "Georgia"        "Hawaii"         "Idaho"         
[13] "Illinois"       "Indiana"        "Iowa"           "Kansas"        
[17] "Kentucky"       "Louisiana"      "Maine"          "Maryland"      
[21] "Massachusetts"  "Michigan"       "Minnesota"      "Mississippi"   
[25] "Missouri"       "Montana"        "Nebraska"       "Nevada"        
[29] "New Hampshire"  "New Jersey"     "New Mexico"     "New York"      
[33] "North Carolina" "North Dakota"   "Ohio"           "Oklahoma"      
[37] "Oregon"         "Pennsylvania"   "Rhode Island"   "South Carolina"
[41] "South Dakota"   "Tennessee"      "Texas"          "Utah"          
[45] "Vermont"        "Virginia"       "Washington"     "West Virginia" 
[49] "Wisconsin"      "Wyoming"

Dados: USArrests

  • Inicialmente, examinamos brevemente os dados. Podemos notar que as variáveis têm médias substancialmente diferentes.
summary(USArrests)
     Murder          Assault         UrbanPop          Rape      
 Min.   : 0.800   Min.   : 45.0   Min.   :32.00   Min.   : 7.30  
 1st Qu.: 4.075   1st Qu.:109.0   1st Qu.:54.50   1st Qu.:15.07  
 Median : 7.250   Median :159.0   Median :66.00   Median :20.10  
 Mean   : 7.788   Mean   :170.8   Mean   :65.54   Mean   :21.23  
 3rd Qu.:11.250   3rd Qu.:249.0   3rd Qu.:77.75   3rd Qu.:26.18  
 Max.   :17.400   Max.   :337.0   Max.   :91.00   Max.   :46.00
  • Vemos que há, em média, três vezes mais estupros (Rape) do que homicídios (Murder), e mais de oito vezes mais agressões (Assault) do que estupros.

Dados: USArrests

  • Também podemos examinar as variância das quatro variáveis usando a função apply():
apply(USArrests, 2, var)
    Murder    Assault   UrbanPop       Rape 
  18.97047 6945.16571  209.51878   87.72916

  • Observe que a função apply() nos permite aplicar uma função — neste caso, a função var() — a cada linha ou coluna da data frame.

  • O segundo argumento aqui indica se desejamos calcular a variância das linhas, \(1\), ou das colunas, \(2\).

  • Não surpreende que a variabilidade seja muito diferente.

Dados: USArrests - Variabilidade

  • A variável UrbanPop mede a porcentagem da população de cada estado que vive em área urbana, que não é um número comparável ao número de estupros em cada estado por 100.000 indivíduos.

  • Se não padronizamos as variáveis antes de realizar a PCA, a maioria dos componentes principais que observamos seriam impulsionados pela variável Assault, uma vez que ela tem de longe as maiores média e variância.

  • Assim, é importante padronizar as variáveis para que tenham média zero e desvio padrão um antes de realizar a PCA.

Análise de Componentes Principais com prcomp()

  • Agora realizamos a PCA usando a função prcomp(), que é uma das várias funções em R que executam PCA:
prcomp_output <- prcomp(USArrests, scale = TRUE)

  • Por padrão, a função prcomp() centraliza as variáveis para terem média zero.

  • Usando a opção scale = TRUE, escalamos o variáveis para que tenham desvio padrão igual a um.

PCA com prcomp()

  • A saída de prcomp() contém uma série de quantidades úteis.
names(prcomp_output)
[1] "sdev"     "rotation" "center"   "scale"    "x"
  • Os componentes center e scalce correspondem às médias e desvios padrão das variáveis que foram padronizadas antes da implementação da PCA.
prcomp_output$center
  Murder  Assault UrbanPop     Rape 
   7.788  170.760   65.540   21.232 
prcomp_output$scale
   Murder   Assault  UrbanPop      Rape 
 4.355510 83.337661 14.474763  9.366385

PCA com prcomp()

  • A matriz rotation contém as cargas dos componentes principais; cada coluna desta matriz contém o correspondente vetor de cargas do componente principal.
prcomp_output$rotation
                PC1        PC2        PC3         PC4
Murder   -0.5358995 -0.4181809  0.3412327  0.64922780
Assault  -0.5831836 -0.1879856  0.2681484 -0.74340748
UrbanPop -0.2781909  0.8728062  0.3780158  0.13387773
Rape     -0.5434321  0.1673186 -0.8177779  0.08902432

PCA com prcomp()

  • Esta matriz é chamada de matriz de rotação, porque quando multiplicamos a matriz \(\bf X\) pela matriz de rotação, obtemos as coordenadas dos dados no sistema de coordenadas rotacionado. Essas coordenadas são os scores dos componentes principais.
prcomp_output$rotation
                PC1        PC2        PC3         PC4
Murder   -0.5358995 -0.4181809  0.3412327  0.64922780
Assault  -0.5831836 -0.1879856  0.2681484 -0.74340748
UrbanPop -0.2781909  0.8728062  0.3780158  0.13387773
Rape     -0.5434321  0.1673186 -0.8177779  0.08902432

PCA com prcomp()

  • Vemos que há 4 componentes principais distintos. Isso era esperado, porque em geral, existem \(\min(n-1,p)\) componentes principais informativos em um conjunto de dados com \(n\) observações e \(p\) variáveis.

  • Usando a função prcomp(), não precisamos multiplicar explicitamente os dados pelos vetores de cargas do componente principal para obter os vetores de scores do componente principal.

  • Em vez disso, a matriz x (\(50 \times 4\)) tem como colunas os vetores de scores dos componentes principais. Ou seja, a \(k\)ésima coluna de x é o \(k\)-ésimo vetor de score do componente principal.

PCA com prcomp()

  • Visualizando as dimensões da matriz de rotação x e os vetores de scores dos componentes principais, que são as colunas de x:
dim(prcomp_output$x)
[1] 50  4
prcomp_output$x
                       PC1         PC2         PC3          PC4
Alabama        -0.97566045 -1.12200121  0.43980366  0.154696581
Alaska         -1.93053788 -1.06242692 -2.01950027 -0.434175454
Arizona        -1.74544285  0.73845954 -0.05423025 -0.826264240
Arkansas        0.13999894 -1.10854226 -0.11342217 -0.180973554
California     -2.49861285  1.52742672 -0.59254100 -0.338559240
Colorado       -1.49934074  0.97762966 -1.08400162  0.001450164
Connecticut     1.34499236  1.07798362  0.63679250 -0.117278736
Delaware       -0.04722981  0.32208890  0.71141032 -0.873113315
Florida        -2.98275967 -0.03883425  0.57103206 -0.095317042
Georgia        -1.62280742 -1.26608838  0.33901818  1.065974459
Hawaii          0.90348448  1.55467609 -0.05027151  0.893733198
Idaho           1.62331903 -0.20885253 -0.25719021 -0.494087852
Illinois       -1.36505197  0.67498834  0.67068647 -0.120794916
Indiana         0.50038122  0.15003926 -0.22576277  0.420397595
Iowa            2.23099579  0.10300828 -0.16291036  0.017379470
Kansas          0.78887206  0.26744941 -0.02529648  0.204421034
Kentucky        0.74331256 -0.94880748  0.02808429  0.663817237
Louisiana      -1.54909076 -0.86230011  0.77560598  0.450157791
Maine           2.37274014 -0.37260865  0.06502225 -0.327138529
Maryland       -1.74564663 -0.42335704  0.15566968 -0.553450589
Massachusetts   0.48128007  1.45967706  0.60337172 -0.177793902
Michigan       -2.08725025  0.15383500 -0.38100046  0.101343128
Minnesota       1.67566951  0.62590670 -0.15153200  0.066640316
Mississippi    -0.98647919 -2.36973712  0.73336290  0.213342049
Missouri       -0.68978426  0.26070794 -0.37365033  0.223554811
Montana         1.17353751 -0.53147851 -0.24440796  0.122498555
Nebraska        1.25291625  0.19200440 -0.17380930  0.015733156
Nevada         -2.84550542  0.76780502 -1.15168793  0.311354436
New Hampshire   2.35995585  0.01790055 -0.03648498 -0.032804291
New Jersey     -0.17974128  1.43493745  0.75677041  0.240936580
New Mexico     -1.96012351 -0.14141308 -0.18184598 -0.336121113
New York       -1.66566662  0.81491072  0.63661186 -0.013348844
North Carolina -1.11208808 -2.20561081  0.85489245 -0.944789648
North Dakota    2.96215223 -0.59309738 -0.29824930 -0.251434626
Ohio            0.22369436  0.73477837  0.03082616  0.469152817
Oklahoma        0.30864928  0.28496113  0.01515592  0.010228476
Oregon         -0.05852787  0.53596999 -0.93038718 -0.235390872
Pennsylvania    0.87948680  0.56536050  0.39660218  0.355452378
Rhode Island    0.85509072  1.47698328  1.35617705 -0.607402746
South Carolina -1.30744986 -1.91397297  0.29751723 -0.130145378
South Dakota    1.96779669 -0.81506822 -0.38538073 -0.108470512
Tennessee      -0.98969377 -0.85160534 -0.18619262  0.646302674
Texas          -1.34151838  0.40833518  0.48712332  0.636731051
Utah            0.54503180  1.45671524 -0.29077592 -0.081486749
Vermont         2.77325613 -1.38819435 -0.83280797 -0.143433697
Virginia        0.09536670 -0.19772785 -0.01159482  0.209246429
Washington      0.21472339  0.96037394 -0.61859067 -0.218628161
West Virginia   2.08739306 -1.41052627 -0.10372163  0.130583080
Wisconsin       2.05881199  0.60512507  0.13746933  0.182253407
Wyoming         0.62310061 -0.31778662  0.23824049 -0.164976866

PCA: biplot

  • Podemos representar graficamente os dois primeiros componentes principais com um biplot():
biplot(prcomp_output, scale = 0)

  • O argumento scale = 0 de biplot() garante que as setas sejam dimensionadas para representar as cargas;

  • outros valores para scale fornecem biplots ligeiramente diferentes com interpretações diferentes.

PCA: biplot

  • Observe que este gráfico é uma imagem espelhada da Figura 2 dos slides da aula de 18/10.

  • Lembre-se de que os componentes principais são únicos à exceção de uma mudança de sinal, assim, vamos tentar reproduzir a Fig. 2 fazendo algumas pequenas alterações:

prcomp_output$rotation = -prcomp_output$rotation
prcomp_output$x = -prcomp_output$x
biplot(prcomp_output, scale = 0, cex = 0.5)

PCA: biplot

PCA: Porcentagem da Variância Explicada

  • A função prcomp() também gera o desvio padrão de cada componente principal.

  • Por exemplo, para os dados USArrests, podemos acessar esses desvios padrão com:

prcomp_output$sdev
[1] 1.5748783 0.9948694 0.5971291 0.4164494
  • A variância explicada por cada componente principal é obtida elevando os desvios padrão ao quadrado:
cp.var <- prcomp_output$sdev^2
cp.var
[1] 2.4802416 0.9897652 0.3565632 0.1734301

PCA: Porcentagem da Variância Explicada

  • Para calcular a proporção da variância explicada por cada componente principal, simplesmente dividimos a variância explicada por cada componente pela variância total explicada por todos os quatro componentes principais:
pve <- cp.var / sum(cp.var)
pve
[1] 0.62006039 0.24744129 0.08914080 0.04335752
  • Vemos que o primeiro componente principal explica \(62,0\,\%\) da variação nos dados, o segundo componente principal explica \(24,7\,\%\) da variação, e assim por diante.

PCA: Porcentagem da Variância Explicada

  • Podemos plotar a PVE explicada por cada componente, bem como a PVE acumulada, com:
par(mfrow = c(1, 2))

plot(pve, 
     xlab = "Componente Principal",
     ylab = "Proporção oda Variância Explicada", 
     ylim = c(0, 1),
     type = "b",
     col = "blue")

plot(cumsum(pve), 
     xlab = "Componente Principal",
     ylab = "Proporção Acumulada da Variância Explicada",
     ylim = c(0, 1), 
     type = "b",
     col = "blue")

PCA: Porcentagem da Variância Explicada

PCA: Porcentagem da Variância Explicada

  • Mas, podemos obter os devios padrão e a proporção da variância explicada por cada componente aplicando a função summary ao resultado de prcomp:
prcomp_output <- prcomp(USArrests, scale = TRUE)
summary(prcomp_output)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.5749 0.9949 0.59713 0.41645
Proportion of Variance 0.6201 0.2474 0.08914 0.04336
Cumulative Proportion  0.6201 0.8675 0.95664 1.00000

Package: FactoMineR

  • Pacote que disponibiliza métodos exploratórios de análise de dados para resumir, visualizar e descrever dados multivariados.

  • Os principais métodos de PCA são disponibilizados:

    • Análise de Componentes Principais (PCA) quando as variáveis são quantitativas,
    • Análise de Correspondência (CA) e Análise de Correspondência Múltipla
      (MCA) quando as variáveis são categóricas,
    • Análise de Fatores Múltiplos quando as variáveis são estruturadas em grupos, etc. e;
    • Aálise de Clusters: método hierárquico

Package: facotextra

  • Fornece algumas funções simpes de usar para extrair e visualizar o resultado de análises de dados multivariados.

  • Compatível com os resultados do pacote FactoMineR

PCA usando o pacote FactoMineR

  • Argumentos da funcão FactoMineR::PCA:

  • scale.unit = TRUE: por padrão as variáveis são centradas (\(x_i - \bar{x}\)), se TRUE, divide-se as variáveis centradas pelo desvio-padrão: \(\frac{(x_i - \bar{x})}{s_{x_i}}\).

  • ncp define o número máximo de componentes: \(\min(n-1,p)\)

PCA usando o pacote FactoMineR 

library(FactoMineR)

fmine_output = FactoMineR::PCA(USArrests, scale.unit = TRUE, ncp = 4, graph = F)
summary(fmine_output)

Call:
FactoMineR::PCA(X = USArrests, scale.unit = TRUE, ncp = 4, graph = F) 


Eigenvalues
                       Dim.1   Dim.2   Dim.3   Dim.4
Variance               2.480   0.990   0.357   0.173
% of var.             62.006  24.744   8.914   4.336
Cumulative % of var.  62.006  86.750  95.664 100.000

Individuals (the 10 first)
                Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
Alabama     |  1.574 |  0.986  0.783  0.392 | -1.133  2.596  0.518 |  0.444
Alaska      |  3.051 |  1.950  3.067  0.409 | -1.073  2.327  0.124 | -2.040
Arizona     |  2.089 |  1.763  2.507  0.712 |  0.746  1.124  0.127 | -0.055
Arkansas    |  1.149 | -0.141  0.016  0.015 | -1.120  2.534  0.950 | -0.115
California  |  3.037 |  2.524  5.137  0.690 |  1.543  4.811  0.258 | -0.599
Colorado    |  2.114 |  1.515  1.850  0.513 |  0.988  1.971  0.218 | -1.095
Connecticut |  1.860 | -1.359  1.489  0.534 |  1.089  2.396  0.343 |  0.643
Delaware    |  1.184 |  0.048  0.002  0.002 |  0.325  0.214  0.075 |  0.719
Florida     |  3.070 |  3.013  7.321  0.964 | -0.039  0.003  0.000 |  0.577
Georgia     |  2.366 |  1.639  2.167  0.480 | -1.279  3.305  0.292 |  0.342
               ctr   cos2  
Alabama      1.107  0.080 |
Alaska      23.343  0.447 |
Arizona      0.017  0.001 |
Arkansas     0.074  0.010 |
California   2.010  0.039 |
Colorado     6.726  0.268 |
Connecticut  2.321  0.120 |
Delaware     2.897  0.368 |
Florida      1.866  0.035 |
Georgia      0.658  0.021 |

Variables
               Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
Murder      |  0.844 28.719  0.712 | -0.416 17.488  0.173 |  0.204 11.644
Assault     |  0.918 34.010  0.844 | -0.187  3.534  0.035 |  0.160  7.190
UrbanPop    |  0.438  7.739  0.192 |  0.868 76.179  0.754 |  0.226 14.290
Rape        |  0.856 29.532  0.732 |  0.166  2.800  0.028 | -0.488 66.876
              cos2  
Murder       0.042 |
Assault      0.026 |
UrbanPop     0.051 |
Rape         0.238 |

FactoMineR: Scree Plot

library(factoextra)

fviz_screeplot(fmine_output, choice = "variance", addlabels = TRUE, 
               ylim = c(0, 100))

Package: FactoMineR - Gráficos

  • A função plot.PCA produz dois gráficos:

    • o mapa de fatores das variáveis (colunas).
    • o mapa de fatores das indivíduos (observações), e;

  • Mapa de fatores das variáveis:

    • O mapa de fatores das variáveis apresenta uma visão da projeção das variáveis observadas no plano abrangendo os dois primeiros componentes principais.

    • Isso nos mostra a relação estrutural entre as variáveis e os componentes, e nos ajuda a nomear os componentes.

    • A projeção de um vetor das variável no eixo do componente nos permite ver diretamente a correlação entre a variável e o componente.

    • A ideia desse gráfico é mostrar com qual direção (componete) as variáveis sào correlacionadas. O eixo que representa Dim 1 e Dim 2 contém o coeficiente de correlação de Pearson (\(-1 \leq r \leq +1\)).

FactoMineR: Mapa de fatores das variáveis

plot.PCA(fmine_output, axes=c(1, 2), choix="var")

factoextra: Mapa de fatores das variáveis

fviz_pca_var(fmine_output, col.var = "orange")

FactoMineR: Mapa de fatores das variáveis

  • O mapa de fatores dos indivíduos exibe os indivíduos prjetados sobre os scores (z) dos componentes principais:
plot.PCA(fmine_output, axes = c(1, 2), choix = "ind")

factoextra: Mapa de fatores das variáveis e indivíduos

fviz_pca_biplot(fmine_output, 
                label = "all", 
                col.var = "orange",
                col.ind = "blue",,
                repel = TRUE)

factoextra: Mapa de fatores das variáveis e indivíduos